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Abstract 

Esophageal transport is a physiological process that mechanically transports an 
ingested food bolus from the pharynx to the stomach via the esophagus, a multi¬ 
layered muscular tube. This process involves interactions between the bolus, the 
esophagus, and the neurally coordinated activation of the esophageal muscles. In this 
work, we use an immersed boundary (IB) approach to simulate peristaltic transport 



in the esophagus. The bolus is treated as a viscous fluid that is actively transported 
by the muscular esophagus, which is modeled as an actively contracting, fiber- 
reinforced tube. A simplified version of our model is verified by comparison to an 
analytic solution to the tube dilation problem. Three different complex models of the 
multi-layered esophagus, which differ in their activation patterns and the layouts of 
the mucosal layers, are then extensively tested. To our knowledge, these simulations 
are the first of their kind to incorporate the bolus, the multi-layered esophagus 
tube, and muscle activation into an integrated model. Consistent with experimental 
observations, our simulations capture the pressure peak generated by the muscle 
activation pulse that travels along the bolus tail. These fully resolved simulations 
provide new insights into roles of the mucosal layers during bolus transport. In 
addition, the information on pressure and the kinematics of the esophageal wall due 
to the coordination of muscle activation is provided, which may help relate clinical 
data from manometry and ultrasound images to the underlying esophageal motor 
function. 

Key words: fluid-structure interaction, immersed boundary method, esophageal 
transport, muscle activation 


1 Introduction 


Interactions between fluids and deformable structures are widespread in bio¬ 
logical systems, and such interactions often involve complex moving interfaces 
and large structural deformations [1]. Esophageal transport is one such pro¬ 
cess, whereby the food bolus is transported to the stomach via the esophagus. 
The esophagus is a flexible, multi-layered tube that consists of mucosal, in¬ 
terfacial, circumferential, and longitudinal muscle layers. The pumping force 
required to produce this peristaltic transport process is generated by neurally 
coordinated muscle activation along the esophagus [2,3], and accounting for 
the full physiological details of the transport process is challenging. Simplified 
analytical models can provide some insights into this biophysical process [4] 
but are often limited in their scope. More complete models are needed to in¬ 
vestigate esophageal pathophysiology, such as motility disorders, and hold the 
potential to advance diagnoses and patient treatment. 

Current studies on the modeling of esophageal transport have focused on 
specific, albeit important, subproblems, such as characterizing the material 
properties of each layer of the esophagus tube [5,6,7,8,9,10], investigating the 
flow of bolus with specified time-dependent bolus geometry or known lumen 
pressure [4,11,12], and estimating the muscle active tension based on known 
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time-dependent pressure distribution [13,14], To the best of our knowledge, 
however, there is presently no computational model of esophageal transport 
that couples models of the bolus, the esophageal structure, and the muscle 
activation within an integrative numerical model. 

This work presents one such integrative model of esophagael transport that is 
based on the immersed boundary (IB) method [15], The IB method is an ap¬ 
proach to modeling fluid-structure interaction that was introduced to simulate 
the fluid dynamics of heart valves [16,17], and which has subsequently been 
applied to a broad range of problems in biology [15], The IB method uses an 
Eulerian description of the momentum and incompressibility of the coupled 
fluid-structure system along with a Lagrangian description of the structural 
forces produced by the elasticity or active tension generation of the structure. 
The primary advantage of this formulation is that it avoids the need to em¬ 
ploy body-fitted grids, and thereby eliminates the need to develop complex 
remeshing strategies as the immersed structure deforms [15,18]. In this paper, 
we present a fully resolved fluid-structure intreaction model of esophageal 
transport, which includes detailed descriptions of the esophageal wall, the bo¬ 
lus, and their interaction. The model is fully resolved in the sense that it does 
not assume simplified fluid dynamics or structural deformations. 

In our model, the majority of the computational domain is occupied by the 
immersed body, with a fluid (i.e., the bolus) confined in a narrow lumen. 
Thus, it is important to describe the mechanical response of the esophagus 
tube. To that end, we discretize the continuous fibers of the esophagus into 
springs and beams and associate a volumetric patch with each such spring 
and beam. This allows us to compute the spring and beam parameters from 
the material properties of the esophagus (e.g., its Young’s modulus). To test 
the fiber-based esophagus tube model, we simulate the problem of dilation of 
a three-dimensional tube and compare the numerically obtained inner fluid 
pressure with an analytically derived solution. 

To simulate bolus transport in a physiologically realistic manner, an esophagus 
model is constructed as a four-layered structure. Muscle activation that results 
in the peristaltic motion of the esophagus is modeled via springs with dynamic 
rest lengths. To handle the numerical challenge arising from large deformations 
of the mucosal layer, we employ a locally refined structural discretization that 
ensures that the Lagrangian structure does not “leak” even under very large 
deformations [15]. We consider three cases with different muscle activation 
models and mucosal layer fiber arrangements, and we discuss the key features 
related to muscle cross-section area and pressure peaks during bolus transport. 
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2 Mathematical formulation 


2.1 The immersed boundary method 


The IB formulation of problems of fluid-solid interaction employs an Eulerian 
description for the momentum equation and the divergence-free condition and 
a Lagrangian description of the deformation of the immersed structure and 
the resulting structural forces. Here, we use the same notation for the Eule¬ 
rian and Lagrangian coordinates as detailed in Griflith [19]. Specifically, we 
let X = (a:i,X2,X3) C ri denote fixed Cartesian coordinates, in which Q C 
denotes the fixed domain occupied by the entire fluid-structure system. We 
use s = (si, S 2 , S 3 ) C U to denote the Lagrangian coordinates attached to the 
immersed structure, in which C denotes the fixed material coordinate 
system attached to the structure. Eor simplicity of implementation, we con¬ 
sider that the fluid-structure system possesses a uniform mass density p and 
dynamic viscosity p. This simplification implies that the immersed structure 
is neutrally buoyant and viscoelastic rather than purely elastic. An extension 
of the present mathematical formulation to problems with nonuniform mass 
densities or viscosities is also feasible; see Ref. [20] for details. 

The equations of motion of the coupled fluid-structure system are [15] 

= — Vp(x, f) -I- /iV^u(x, t) 

+ g(x, t), 

V • u(x, t) = 0, 

g(x,f)= [ G(s,t)(J(x - X(s,t))ds, 

Jci 

5X c 

-^(s, t) = u(x, t) S{x - X(s, t)) dx. 

G{s,t) = g[X{s,t)]. 


P f^(x, t) + u(x, t) ■ Vu(x, t) 


( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


Eqs. (1) and (2) are the incompressible Navier-Stokes equations written in 
the Eulerian form, u(x, t) is the Eulerian velocity, p(x, t) is the pressure, and 
g(x, t) is the Eulerian elastic force density. Eq. (5) describes the elastic force 
in the immersed body in Lagrangian form, in which G(s,t) is the elastic 
force density and ^ : X i-)- G is a time-dependent functional that determines 
the Lagrangian force density from the current configuration of the immersed 
structure. Interactions between Lagrangian and Eulerian variables in eqs. (3) 
and (4) are mediated by integral transforms with a three-dimensional Dirac 
delta function kernel 5(x) = Specifically, eq. (3) converts the La¬ 

grangian force density G(s, t) into an equivalent Eulerian force density g(x, t), 
and eq. (4) determines the physical velocity of each Lagrangian material point 
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from the Eulerian velocity field, thereby effectively imposing the no-slip condi¬ 
tion along the fluid-solid interface. The discretized versions of these equations 
used in this work employ a regularized version of the delta function, denoted 
5/i(x) = for details on the construction of such regularized delta 

functions, see Ref. [15]. The details on the spatial discretization of Eulerian 
fluid system (i.e., eqs. (1) and (2)), Lagrangian-Eulerian interaction equations, 
(i.e., eqs. (3) and (4)), and the temporal discretization of the system of equa¬ 
tions can be found in refs. [19,21]. In the following section, we discuss the 
Lagrangian discretization of eq. (5) to characterize the material elasticity of 
the immersed structure. 


2.2 Material elasticity 


The specific form of the mapping function Q : X g is dictated by the 
model of material elasticity of the immersed body. To that end, Chadwick [22], 
Ohayon and Chadwick [23], and Tozeren [24] proposed the “fluid-fiber” and 
the “fiuid-fiber-collagen” models to characterize the material elasticity of bio¬ 
logical tissues. They assumed the tissues to be an aggregation of elastic fibers 
that are embedded in a soft matrix. These models were used to describe the 
esophageal wall by Nicosia and Brasseur [14]. They considered the muscle layer 
as a family of fibers and discarded the elasticity of the soft matrix. Inspired by 
the success of their model [14], we also ignore the elasticity of the soft ground 
matrix in the esophageal wall, and model all esophageal layers as families of 
continuous fibers embedded in the background fluid. Thus, the material elas¬ 
ticity of the esophagus tube is essentially represented by the fiber elasticity, 
which is approximated using springs and beams in our discretized IB scheme. 
By contrast, Ghosh et al. [13] have modeled the esophageal muscle layer as a 
family of incompressible continuous fibers embedded in a soft isotropic matrix. 
Such models will be explored in our future work. 

We describe the fiber-based material elasticity in terms of a strain-energy 
functional E = £'[X(-,t)]. The corresponding Lagrangian elastic force density 
can be derived by taking the Erechet derivative of E as 

pE[X(-, t)] = -J^ G(s, t) • pX(s, t) ds, (6) 

in which p denotes the perturbation of a quantity. Since the elastic properties 
of the structure are described in terms of a family of elastic fibers that resist ex¬ 
tension, compression and bending, the strain-energy functional E = £'[X(-, t)] 
can be decomposed into a stretching part Eg that accounts for the extension 
and compression of the fibers, and a bending part E^ which accounts for the 
resistance of the fibers to bending, i.e., E = Eg + Ef,. 
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2.3 Elastic springs and beams 


We use a 3D spring network to represent the elastic stretching energy Eg of 
the esophagus tube. Three families of springs are employed for the radial, 
circumferential, and axial fibers of the tube. On the contrary, we only include 
axial beams in esophagus model to account for the bending energy E^, as 
the significant curvature change for the long esophageal tube occurs primarily 
along the axial direction. 

In our implementation of the IB method [19], instead of computing the elastic 
force density, the total elastic force is computed for a Lagrangian node, which is 
then spread to the Eulerian grid (via eq. (3)) to obtain the equivalent Eulerian 
force density. To relate the spring/beam constant with the material property 
of the esophagus (such as its Young’s modulus), we associate a volumetric 
patch with each spring and beam. To obtain the volumetric patch, we describe 
the tube in cylindrical coordinates (r, 6, z) with the reference configuration 
given by (a < r < 6,0 < 6* < 27r,0 < z < 1), where a, b and I are the 
inner radius, the outer radius and the length of the tube in the reference 
configuration, respectively. Let {ri,0j,Zk) denote the node point of the tube, 
with Ar, = rj_|_i — Tj, AOj = Oj^i — Oj and Azk = z^+i—Zk denoting the spacing 
along r, 6 and .2 coordinates, respectively. Then, a case of uniform spacing 
Ar, A6 and Az will be obtained if we let Ar, = Ar, A9j = A9 and Azk = 
Az for the discretization. Depending upon whether any of the three types of 
springs or an axial beam lies in the interior of the structure or its boundary, 
different volumetric patches (interior or boundary patches) will be associated 
with them. We layout the nodal points in such a way that volumetric patches 
of each type Ptype (i-e., patches associated with radial/axial/circumferential 
springs or axial beams) do not overlap with each other and they sum up to 
the volume of the esophagus tube. This can be written as 

f-^type n-Ptype ~ ^ 

(U/ -^type ~ S/ -^type ~ f^sophagus I ^ ^type-i 

in which, Wsophagus is the volume of the esophagus and Wype is the number of 
patches of the same type. 

Eig. 1(a) and 1(b) show the volumetric patches of various springs in (r, d) 
plane and (r, z) plane, respectively. Eig. 2(a) shows the patches associated 
with axial beams. Cases with nonuniform A9 are also considered in this work. 
In particular, we consider cases with A9 for inner layers and 0.5Ad for outer 
layers of the tube. The patches associated with the circumferential springs in 
(r, 9) plane for the nonuniform A9 are shown in Eig. 2(b). Once the patches for 
the discrete springs and beams are obtained, the elastic force due to extension, 
compression and bending of the fibers can be computed with given material 
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parameters. Next we discuss how spring/beam constants are obtained in our 
model from the elastic modulus of the fibers. 


Spring constant: We assume a linear relationship between the stress and strain 
of a spring, sp, which provides resistance to extension and compression of a 
fiber. If the stress and strain of the spring are defined with respect to the 
undeformed configuration, then we have 

a = — = be, 
l-L 


in which a, e, and F are the stress, strain, and internal force of the spring, 
respectively, S is the elastic modulus of the fiber, I is the current length, L is 
the rest length and A is the undeformed cross sectional area. Consequently, 
the expression for the spring nodal forces and for nodes Ii and I 2 
connected by the spring are given by 


( 8 ) 

( 9 ) 


■pP ^ _-ph ^ 
sp sp 


Lj 


X 


h 

sp 


^sp ^sp 


v /2 


X 


h 

sp 


^sp ^sp 


( 10 ) 


in which, K = SA/L is the spring constant (also called the spring stiffness). 
The undeformed cross sectional area A can be obtained from the associated 
patch of the spring, A = Vgp/L, where 14p is the volume of patch of the spring. 


Beam constant: A beam b associated with three nodes p, I 2 , and I 3 , provides 
resistance to bending and sets a preferred (possibly a time-dependent) cur¬ 
vature of the fiber. The nodal forces for the beam F^C F^^, and F^® can be 
obtained from the Frechet derivative of the bending energy as 




ds'^ 


ds"^ 


2 

ds. 


( 11 ) 


in which, Cb is the bending stiffness and P/ is the associated patch of the beam. 

is the preferred configuration, which in present work satisfies = 0. 
Thus, 

Gb ■ pXfc ds = - [ ^ Cb^^- ■ ds. (12) 


'Pi 


'Pi 


ds"^ ds"^ 


We approximate the second derivative (with respect to the arc length s) in 
eq. (12) in the reference configuration of the tube, where the three beam nodes 
are assumed to be on the same line (i.e., have zero curvature). This basically 
implies small bending deformation. With g, (for i = 1,2,3) labeling the three 
nodes, we construct shape function W for node g* in the local coordinate 
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Ni{q) = 

N 2 {q) = 

Nsiq) = 

Thus, eq. (12) evaluates as 


: E|=i Xj^W,(g) and pX, = ELi 


{q - q 2 ){q - q?i) 

d^Ni 2 

(13) 

(a -q 2 ){qi -qs)’ 

dq^ (qi - q 2 ){qi - qs) 

{q-qi){q-q3) 

d‘^N2 2 

(14) 

{q 2 - qi){q 2 - qs)’ 

dq^ (g2 - gi)(g2 - gs) 

{q-q3){q-qi) 

d'^Ns 2 

(15) 

{q3 - qi){q3 - q2)' 

dq^ {qs - qi){q3 - q2) 


iPl 


Gfc • p^b ~ ~ 


d^pXb 


' *^9 0 9 ^q 

Pi oq^ dq‘^ 

— -Vb Cb^b 


= f;* • pxj*. 


dq^ dq^ 


(16) 


in which, is the volume of the patch associated with the beam and = 

. For the case of uniform distance between neighboring 
nodes of a beam in the reference configuration, i.e., when q 2 —qi = qz—q 2 = Ag, 
simple expression for the nodal forces can be obtained as below, 






-Vjcb 

Aq4 


1 



1 -2 1 




xl- 



(17) 


The bending coefficient in the above equation can be determined from the 
relation = ^, in which / and A are the second moment of area and area of 
cross section of the beam, respectively. 


3 Verification case: A 3D dilation problem 

3.1 Problem formulation 


For the fiber-based IB scheme, Griffith [21,25] has studied convergence proper¬ 
ties for various 2D cases. Here we present a test case of a 3D elastic cylindrical 
tube undergoing dilation to verify the solution methodology. A cylindrical tube 
composed of three families of continuous fibers (i.e., axial, circumferential and 
radial fibers) is immersed in a rectangnlar fluid domain. We nondimensional- 
ize the system based on the inner radius of the cylindrical tube, the density 
and viscosity of the fluid. The tube’s reference configuration is described in 
cylindrical coordinates (r, 6 , z) with l<r<l-|-T, 0 < 9 < 2%, 0 < .2 < 10, 
in which T is the tube thickness. The fluid domain is described in Cartesian 



coordinates {x,y,z) with —2<x<2, —2<y<2,0<z<10. The dimen¬ 
sionless Young’s modulus S of all fibers is taken to be the same, 5” = 4 x 10^. 

The dilation process includes two phases. The first phase is an inflation phase. 
This is modeled by adding fluid in the domain from the top end while keeping 
its bottom end closed. The second phase is the relaxation phase. The relax¬ 
ation phase continues until a stationary state is reached, at which the inertial 
and viscous terms are approximately three orders of magnitude smaller than 
the pressure term. Thus, the pressure force from the constraint of fluid incom¬ 
pressibility is balanced by the elastic force in the tube. The simulations are 
carried by specifying the following boundary conditions for the fluid domain: 
traction (normal and tangential) free boundary conditions for the four lateral 
surfaces of the domain; zero-velocity boundary condition for the bottom sur¬ 
face; and a time-dependent velocity boundary condition for the top surface. 
On the inflow surface, the tangential velocity is zero and the normal velocity 
is u-n = —UQf(t)(l — x'^ — y'^) if x'^ + y^ < 1, and zero otherwise. Here, f{t) is a 
decreasing function of time which vanishes at the end of the inflation process. 

We assume a plain strain state for this relatively long tube at the final equi¬ 
librium state. Let U (r) denote the radial displacement field of the tube in the 
middle region, then by assuming plane strain conditions (i.e., a sufficiently 
long tube) in the final equilibrium state, an analytic expression for the inner 
pressure can be obtained (see Appendix A), 


P- — ^ 

inner ^ 


log -h 


rf-C 

r\ 


-S 


log Wo + - C -h 


ri-C 



(18) 


Here, C — rf — Rf, r, and Tq are the deformed inner and outer radius, re¬ 
spectively; Ri and Rq denote the initial inner and outer radius, respectively. 
Thus, Vi = Ri + U (n), with U (n) denoting the radial displacement of the inner 
surface of the tube and Tq = R‘^ — Rf + rf. We use the observed 17(r,) in 

our simulations, and compare the numerical and predicted analytic values of 
the inner pressure, respectively Pnumericai and Panaiytic- 


3.2 Verification results 


Here, we conduct test cases with different tube thickness and dilation levels, 
where higher dilation level is simulated by adding more fluid into the tube 
during the dilation phase. For each case, the Eulerian computational domain 
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is discretized using an AA x x 50 Cartesian grid, whereas the Lagrangian 
structural domain is described using a cylindrical coordinate system (r, z) 
and discretized using an x x 100 mesh. To prevent the fluid from leaking 
out of the structure, we keep the grid size of Lagrangian mesh smaller than 
that of Eulerian mesh. The grid number for cases with different tube thickness 
is listed in Table 1. For a thicker tube of thickness T = 1, we use nonuniform 
A0 for inner and outer layers, with = 64 for the inner and = 128 for 
the outer layers. 

Table 1 

Error in the inner pressure for different tube thickness T. U(ri) and ^numerical 
measured radial displacement of the inner surface and inner pressure in the middle 
section of the tube, respectively. The relative error 6^ = Tnumemcai -Panaiytid 

_’ ^^^ _Cknalytic_ 


T 

Nr X Ne 

N 

U{rO 

-^numerical 

-^analytic 


0.2 

4 X 128 

40 

0.07989 

399.4 

403.14 

9.28e-3 

0.5 

10 X 128 

40 

0.07410 

699.9 

710.78 

1.53e-2 

1 

6 X 64 + 4 X 128 

20 

0.07313 

961.4 

972.30 

1.12e-2 


Table 2 

Error in the inner pressure for different dilation levels with tube thickness T = 0.2. 
Higher dilation level is simulated by adding more fluid in the tube, as shown by the 
increase of U{n). The relative error Cp = I ^numerical-^analytic I 

■^ctricily t ic 


Dilation level 

U{ri) 

-^numerical 

-^analytic 

^p 

Level 1 

0.02827 

159.7 

161.24 

9.57e-3 

Level 2 

0.07989 

399.4 

403.14 

9.28e-3 

Level 3 

0.15658 

652.7 

663.85 

1.68e-2 

Level 4 

0.23483 

805.2 

840.9 

4.25e-2 


As can be seen from Table 1 and Table 2, our fiber-based tube model is able 
to capture the analytical trend for various tube thickness and dilation levels, 
with relative error below 5%. 


4 Esophageal transport 


In the previous section, we showed that our IB formulation for an elastic tube 
(with an appropriate arrangement of fibers) is able to capture the analytical 
trend of a dilation process. In this section, we extend the elastic tube model 
to describe esophageal transport. 
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The human esophagus is a long multi-layered composite tube that consists of 
inner mucosal-submucosal layers (collectively referred to as “mucosal” layer) 
and outer muscle layers which in turn include the circular and longitudinal 
muscle layers (so named because of their hber orientations [26]). The anatomy 
of the human esophagus is illustrated in Fig. 3. In-vitro tests show that there 
exists a weak connecting tissue, referred to as the “interfacial” layer, between 
the muscle and mucosal layers [27]. The bolus (depending on its content) is 
generally considered as a Newtonian fluid, with its viscosity varying from one 
centipoise (cP) to several hundred centipoise [28]. The overall volume of the 
bolus is on the order of a few milliliters from the clinical study [29]. Studies of 
Pouderoux et ah [2] and Mittal et ah [3] suggest that during bolus transport, 
the circular muscle contraction is well coordinated with the longitudinal mus¬ 
cle shortening. It is this key activation pattern that is used in our esophageal 
muscle model that enables the transport of the “tear” shaped bolus [4,11] 
through the esophagus. 


4-1 Geometry, boundary conditions and material properties 


The reference conflguration of the esophagus model is taken to be a long 
straight cylindrical tube made up of elastic libers. There are five important 
components in our esophagus model: (1) inner mucosal (IM) layer; (2) outer 
mucosal (OM) layer; (3) interfacial (IF) layer; (4) circular muscle (CM); and 
(5) longitudinal muscle (LM). The IM and OM layers together represent the 
mucosal layer of the esophagus, which we split into two layers for numerical 
purposes (see Sec. 4.3). The length of the esophagus tube is taken to be 240 
mm, as the typical human esophagus length is in the range of 180-250 mm [30]. 
The thin liquid layer confined in the narrow esophageal lumen is assumed to 
have a circular cross section in the reference configuration with a radius of 0.3 
mm. The thickness of each esophageal wall component is obtained based on the 
clinical data of human esophagus at non-rest state (i.e., with intruded catheter 
in the esophagus). The thickness of each layer at rest is listed in Table 3, based 
on the clinical data of Mittal et ah [31]. The entire esophagus is immersed in a 
fluid region of size (—7 mm, 7 mm) x (—7 mm, 7 mm) x (—25 mm, 245 mm). On 
the six surfaces of the fluid box, we impose stress-free boundary conditions. We 
also fix the esophageal top end, which, in physiological situation, is constrained 
by the upper esophageal sphincter. The overall schematic is shown in Fig. 4. 
We here consider the transport of an initially filled bolus in the upper end of 
esophagus. For all the cases presented here, we take uniform viscosity of 10 
cP and uniform density of Ig/cm for both the fluid and the esophagus tube. 

The esophageal tissue is generally modeled as a nonlinear anisotropic elastic 
or pseudo-elastic material. The reported material properties (such as the mod¬ 
ulus of the muscle layers in the circumferential and longitudinal orientations) 
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are, however, substantially different [5,6,7,8,9,10]. Here we assume an elastic 
behavior of each fibrous esophageal layer. Such models have been used before 
to estimate the muscle tension [14]. The elastic modulus for the radial, cir¬ 
cumferential, and axial fibers is taken as S' = 4 kPa. The interfacial layer that 
loosely connects the muscle and mucosal layer consists of only radial fibers 
with modulus S = 0.0004 kPa. For the mucosal layer, Stavropoulou [10] char¬ 
acterized the elasticity of the mucosal layer in their axial and circumferential 
directions, which can be represented by axially-circumferentially-radially ar¬ 
ranged fiber network. By contrast, Natali et al. [6] considered the mucosal 
layer to comprise two families of helical fibers embedded in a matrix tissue. 
Thus, to demonstrate the capabilities of our modeling approach, we consider 
both fiber arrangements of the mucosal layer in Sections 5.1 and 5.3. 

Table 3 

Thickness of each esophageal layer: clinical data [31] and data used in our model. 
Note that the clinical test obtains the thickness of each layer at its non-rest state, 
with esophageal lumen dilated by the intruded catheter. Computer model adopts the 
thickness of each layer at its rest state. The thickness of mucosal layer measured 
based on the ultrasound image is much lower than the thickness at rest, as the 
intruded catheter will distend the mucosal layer significantly and reduce the layer 
thickness._ 


Unit (mm) 

Lumen 

Mucosa 

IF 

CM 

LM 

radius 

thickness 

thickness 

thickness 

thickness 

Clinical data 

3.5 

1.85 

NA 

0.55 

0.49 

Computer model 

0.3 

3.2 

0.6 

0.6 

0.6 


4-2 Muscle activation 


In-vivo experiments [2,3] show that during normal esophageal transport, there 
is well-coordinated circular muscle (CM) contraction and longitudinal muscle 
(LM) shortening. A quantitative model to characterize the contraction and 
shortening process in terms of neuronal firing or reaction kinetics in muscles 
is not available; however, experiments show that there is a precise synchrony 
between the two types of muscle activation patterns [2,3]. In our model, this 
sequential activation is implemented by dynamically changing the rest lengths 
of springs. Specifically, let 2 : denote the vertical coordinate based on the initial 
configuration of the esophageal tube, with the bottom end of the esophagus 
as the origin z = 0, and the top as the end z = L. Then an active spring 
representing a section of one active muscle fiber has its rest length r{z,t) 


12 



given by 


(ro iit-to < ^ 

r{z, t) = I (1 - a{z, t))ro ii ^ < t - to < ^ ^ (19) 

[ro iit-to> ^ + ^ 


in which, ro is the spring’s initial rest length, c is the speed of the activation 
wave, to is the initiation time of activation, a{z, t) is the reduction ratio, and 
AL is the contracting segment’s length in the reference coordinate system. 
Eq. (19) gives the rest length of a spring at its rest, activation and relax¬ 
ation state, respectively. The equation also shows, at any time, the whole 
esophageal tube has a contracting segment with a vertical length AL. The 
variation of muscle activation along this contracting segment will likely in¬ 
fluence bolus transport. To understand this influence, we propose two muscle 
activation models, namely uniform muscle activation and nonuniform muscle 
activation. They differ by how the reduction ratio a{z., t) is distributed along 
the contracting segment AL: 

Uniform muscle activation: 

a{z,t) = do (20) 

Nonuniform muscle activation: 

a(;2,f) (21) 

where uq is a constant, ZQ{t) is the ^-coordinate at the vertical center of the 
contraction segment, and A is the parameter that controls the width of the 
Gaussian distribution in eq. (21). The common parameters of muscle activation 
model used in all the cases of esophageal transport are listed in Table 4. 
a{z,t) differs in test cases depending on whether uniform muscle activation or 
nonuniform muscle activation is used. 

Table 4 

Model parameters for the circular muscle (CM) contraction and longitudinal muscle 
(LM) shortening used in all the cases. The muscle activation model is based on 
eq. (19). 


Muscle activation type 

c (mm/s) 

AL (mm) 

to (s) 

CM contraction 

100 

60 

0 

LM shortening 

100 

60 

0 
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4-.3 Numerical issues 


Esophageal transport involves multiple length scales, which is evidenced by 
the fact that the esophageal length is 240 mm, while the lumen radius at rest is 
only 0.3 mm. The requirement of resolving the narrow lumen dictates the grid 
size of the problem. More challenging in terms of computational modeling 
is the large deformation of the inner lumen, due to the dilation caused by 
bolus movement. This can be seen in Fig. 7 (Section 5.1). For good numerical 
resolution of the transport problem, we use Ax = Ay = 0.2 mm and Az = 
1 mm for the Eulerian mesh. The choice of the Lagrangian mesh size must also 
address the issues of large deformations of lumen and the external fluid leaking 
into the tube, which could occur when the Lagrangian mesh becomes coarser 
than its Eulerian counterpart [15]. In our implementation, considering the 
fact that significant deformation is mainly confined in the inner mucosa (IM) 
layer, we employ a refined mesh in the IM layer along the circumferential and 
axial orientations to form an “impermeable” surface. For the outer layers with 
relatively small dilation, a relatively coarser Lagrangian mesh is used to reduce 
the computational cost. The mesh sizes for various esophageal components are 
fisted in Table 5. The axial beams included in the model provide resistance to 
curvature changes in the axial direction that are associated with buckling of 
the tube. The time step At needs to satisfy the stability constraints from both 
the fluid and solid system. Based on empirical tests, we choose At = 0.02 ms. 
The total time for the transport is about 2.5 s, which requires about 125,000 
timesteps. The relative change in the bolus volume is within 0.8% before the 
bolus begins to empty from the bottom of the esophagus. This indicates that 
the immersed esophagus model is relatively “water-tight”. 

Table 5 

Grid size along r, 6 ,z orientations (denoted as Ar,rA 0 , Az) for each layer of the 
esophagus in the reference configuration. Note that the grid size in circumferential 
orientation rAg for each layer should be understood as an average value, since 
the coordinate r increases from the inner surface to the outer surface while is 
constant in each layer. 


Grid size (0.1mm) 

IM 

OM 

IF 

CM 

LM 

A^ 

2 

4 

3 

3 

2 

rAe 

0.5 

2 

1.9 

2.2 

2.5 

A, 

4 

8 

8 

8 

8 
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5 Results 


5.1 Case 1: Axially-circumferentially-radially arranged mueosal fibers with 
uniform muscle aetivation 


Here we report a case study of the esophageal transport that considers the 
mucosal layer to be composed of axially-circumferentially-radially arranged 
fibers. Experiments show the intact mucosal layer is highly folded at rest (see 
Fig. 1 in Ref. [8]), so we take a relatively low moduli for circumferential, 
radial, and axial fibers as 0.004 kPa, 0.004 kPa and 0.04 kPa, respectively, to 
qualitatively capture its effective elastic property. We use uniform activation 
model as described in eq. (20). The reduction ratio oq for CM contraction and 
LM shortening is taken to be 0.6 and 0.5, respectively. Other parameters are 
fisted in Table 4. 

As shown in the Fig. 5, the bolus is transported to and emptied through the 
bottom of the esophagus in about 2.4 seconds. The running bolus (indicated 
by the negative axial velocity) is confined by the inner mucosal layer. This 
underlines the role of mucosal layers in preparing the “tear-drop” shape of 
the bolus. The pressure distribution shown in Fig. 6 implies that the primary 
pumping force behind the bolus is generated by muscle activation, evidenced 
by the peak pressure at the contraction region. This is consistent with re¬ 
ported experimental observations [3]. We remark that the “tear-drop” shape 
of the running bolus is not specified, but rather is a consequence of the flnid- 
structnre-muscle activation interaction. This is different from previons models 
for bolus transport [4,12,14] where the bolus shape was pre-defined. 

Detailed information on the deformation of each esophageal layer is illustrated 
in Fig. 7. It can be seen that a typical esophageal segment (at each axial lo¬ 
cation) passes through four distinct stages: (a) the segment is at rest; (b) the 
segment is dilated by the incoming bolus; (c) the segment contracts as a re¬ 
sult of the incoming activation wave; and (d) the segment relaxes after the 
activation wave passes. Figs. 5 and 7 show that the mncosal layer plays an 
important role in the bolus transport. First, the inner layer of mucosa shapes 
the running bolus by dynamically closing (or narrowing) the lumen above the 
bolus region while opening the Inmen from below. Second, the pronounced 
axial movement of mucosal layer “lubricates” the running bolus, which helps 
to achieve better transport efficiency. Previous studies on bolus transport ex¬ 
cluded mucosal layers, and did not capture this important feature of bolus 
movement [4,12,14]. Quantitative results for the cross-sectional area (CSA) 
of the esophageal layers, the bolus, as well as for the pressure distribution 
along the Inmen center at time instant t = 1.2 s are shown in Fig. 8. At the 
contraction region behind the bolus, a pressure peak and an increased muscle 
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CSA (i.e. the sum of CSA of LM layer and CSA of CM layer) coexist, which 
is consistent with the experimental observation of Mittal et al. [3]. 


5.2 Case 2: Axially-circumferentially-radially arranged mueosal fibers with 
non-uniform muscle activation 


Here we present the second case study of esophageal transport. We use the 
same geometry and material model of esophagus as of Case 1 in Section 5.1, 
except that we use nonuniform muscle activation model. We take A as 10 mm 
for both LM shortening and CM contraction, and we set reduction ratio Oq for 
CM contraction and LM shortening as 0.6 and 0.5, respectively. Other param¬ 
eters used in the muscle activation model are listed in Table 4. The transport 
phenomenon is shown in Fig. 9 and Fig. 10. It is clear that the nonuniform 
muscle activation results in more pronounced contraction, as illustrated in 
Fig. 11. We remark that a distinctive muscle CSA peak overlaps with the 
pressure peak in this case, as shown in Fig. 12. This is different from Case 1 in 
Section 5.1, which resulted in a plateau of increased muscle CSA (see Fig. 8). 
A coexisting muscle CSA and pressure peak is also observed in the clinical 
test of Mittal et al. [3], which, in conjunction with our simulation, implies 
that a synchronous nonuniform LM shortening and CM contraction probably 
corresponds to the normal physiology of esophageal muscle activation. 


5.3 Case 3: Helical mucosal fibers with uniform muscle activation 


Natali et al. [6] proposed a constitutive model of a multi-layered esophagus, 
in which the mucosal layer consists of two families of helical fibers, and the 
muscle layer is composed of axially and circumferentially arranged fibers. Here 
we present a case with such a helical mucosal fiber arrangement. Two families 
of helical fibers run in the (r, 0) surfaces, as shown in Fig. 13. Radial springs 
are used to link the (r, 6 ) surfaces. The modulus of the helical fibers is taken 
to be the same as that of fibers in the muscle layer, while the modulus of 
the radial fibers of the mucosal layer are one order of magnitude weaker than 
the modulus of helical fibers. We use uniform muscle activation model with 
reduction ratio oq for CM contraction and LM shortening as 0.5 and 0.3, 
respectively. Other parameters of muscle activation model are listed in Table 4. 
Fig. 14 and Fig. 15 shows the bolus transport for this case. A high pressure 
region exists along the contraction segment, which pushes the bolus down, 
similar to what is observed in Case 1 in Section 5.1. The details of deformation 
are illustrated in Fig. 16, which also shows the four distinct stages. However, 
the deformation pattern is more regular than the previous cases, which may 
be attributed to the helical configuration of the mucosal layer. As Fig. 17 
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shows, high pressure and increased muscle CSA overlaps, which indicates a 
synchrony between CM contraction and LM shortening. However, no distinct 
pressure peak exists, which is different from results of Case 1 in Section 5.1. 
This is probably because the helical mucosal layer is less “squeezed”, which is 
evidenced by the slightly decreased CSA. 


6 Remarks on the fiber-based immersed boundary method 


In the previous three sub-sections, we showed results for esophageal transport 
with new insights including roles of mucosal layer and information on lumen 
pressure and kinematics resulting from the synchrony between CM contraction 
and LM shortening, in spite of the complexity of the configuration. Here, we 
add remarks on the application of the classical fiber-based IB method [15] to 
problems involving significant deformations. First, it is noted that the choice of 
number of Lagrangian points per Eulerian grid is challenging. Two Lagrangian 
grid points per Eulerian grid are recommended in each direction to stop the 
fluid from leaking through the immersed structure [15]. However, the esopha¬ 
gus dilates significantly during the bolus transport. It is important to ensure 
that there are at least two grid points per Eulerian grid in each direction in 
the dilated state. This leads to a practical requirement that there should be 
more than two Lagrangian points per Eulerian grid in each direction in the 
rest and relaxation state. Increasing Lagrangian grid refinement relative to 
the Eulerian grid does not necessarily imply higher accuracy, but instead over 
constrains the system. Hence, we need to verify that the results are reasonable 
(Eig. 8). The fact that the number of Lagrangian grid points are determined 
via simulation trials is an undesirable feature. 

One remedy, not within the scope of this work, would be to use an adap¬ 
tive Eulerian mesh, where the adaptive criterion is based on Lagrangian grid 
spacing. Current adaptive meshing of the Eulerian grid is based on the mag¬ 
nitude of velocity gradients. Another approach could be to have an adaptive 
Lagrangian discretization such that the number of Lagrangian grid points in 
the rest and dilated states would differ if the Eulerian grid size is fixed. 

The second issue is related to spurious deformations at the scale of the La¬ 
grangian grid. Multiple Lagrangian grid points per Eulerian grid can result 
in spurious deformation modes at the Lagrangian grid scale. These modes are 
internal to the Lagrangian grid and thus not resolved by the fluid solution 
which is at the scale of the Eulerian grid. As a result, in case of esophageal 
transport, the relaxed configuration recovered after a bolus has passed, has 
residual deformations on the Lagrangian grid scale. These deformations are 
pronounced in the compliant mucosal layer. 
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The third issue pertains to spurious deformations leading to errors in the 
incompressibility of the Lagrangian structure. The incompressibility constraint 
is imposed on the Eulerian grid scale. This is not strictly sufficient to impose 
incompressibility on the Lagrangian grid which is at a sub-grid scale with 
respect to the Eulerian grid. 

One remedy for both issues two and three is to use a spring conhguration with 
diagonal springs in addition to springs oriented in the three orthogonal direc¬ 
tions or to use a helical fiber configuration. Such an approach has been used 
in prior studies involving swimming of two-dimensional eels [32,33]. Another 
remedy for these issues could be to use finite element based Lagrangian im¬ 
mersed structure instead of a fiber-based structure [34]. This is being pursued 
by us but it is not within the scope of this work. 


7 Conclusions 


In this work, we introduce a method based on the volumetric patch to char¬ 
acterize the elasticity of the immersed fiber-based structure employed in the 
typical IB method [15,19,21]. Model verification is performed via comparisons 
between the computational results and an analytic solution to an idealized 
tube dilation problem. Low relative errors (5%) are obtained. To study the 
esophageal transport, we develop a fully resolved active musculo-mechanical 
model that is able to incorporate the liquid bolus, multi-layered esophageal 
wall, and muscle activation into a unified model. We present three cases of 
the esophageal transport that differ in the choice of muscle activation model 
and mucosal fiber arrangement, thereby demonstrating the capabilities and 
generality of the model. The key feature of bolus transport observed experi¬ 
mentally is a “tear-drop” bolus driven by a muscle contraction wave. This is 
captured in our simulations. Moreover, new insights are also provided by fully 
resolved simulations. The simulations show that perfect synchrony between 
LM shortening and CM contraction leads to an overlap of the high pressure 
region and an increased muscle CSA. This helps to relate clinical test data 
from manometry and ultrasound image to the underlying neurally-controlled 
activity, such as the coordination between CM contraction and LM shortening. 
Detailed information on kinematics elucidates the role of the mucosal layer in 
shaping the bolus. Specifically, the mucosal layer provides distensibility to the 
esophageal lumen, and lubricates the running bolus. We remark that this is 
the first study that directly looks at the interaction between the bolus and the 
mucosal layer. Future work should include a parametric study of the effect of 
changing mucosal property on bolus transport. This will help to understand 
certain esophageal diseases that are related to the inflammation of the mucosa. 
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Appendix A: Analysis of a fiber-based tube dilation problem 


Here we derive eq. (19) for the tube dilation problem presented in Section 3. 
The governing equations for an incompressible fluid are 


= V • (Tf = — Vp + /iV^u (22) 

V • u = 0, (23) 

in which cTf is the fluid stress, u is the fluid velocity, p is the fluid density, 
pL is the fluid viscosity and p is the pressure imposing the incompressibility 
constraint. For a neutrally buoyant incompressible viscoelastic structure the 
governing equations are 



p -I- U • Vu^ = V ■ CTe = — Vpe + V • cr -I- (24) 

V • u = 0, (25) 

in which cr^ is the stress tensor in the elastic structure, Pe is the pressure 
that imposes the incompressibility constraint in the elastic structure, a is the 
unknown deviatoric part of the elasticity tensor. The continuity of traction 
and velocity at the fluid-solid interface implies 


(—pi -I- pVu) • n = (— Pel + cr -I- pVu) • n, (26) 

in which n denotes the outward normal vector to the structure interface. The 
above equation implies that there is generally a pressure discontinuity at the 
fluid-structure interface. For the dilation problem at steady state, eq. (22) 
implies a uniform pressure in the tube, denoted as Pmner, since the inertial 
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and viscous terms vanish. Thus, eq. (26) becomes 


Pinner = H • (pi) • H = 11 • (pgl “ Ct) • 11. ( 27 ) 

For the hber-based structure, let R and r denote the initial and deformed 
radial coordinates, respectively, at each material point. At the boundaries of 
the tube, let Pi and Pq denote the initial inner and outer radius, and and 
To denote the current inner and outer radius, respectively. Then we can obtain 
the relationship: P = P(r), based on volume conservation and plain-strain 
deformation: 


P(r) = \Jr‘^ — rf + Rf (28) 

The elasticity of the hber-based tube is represented by three families of springs: 
circumferential, radial, and axial springs. In the (r, 0) plane, only deviatoric 
stresses along the radial and circumferential orientations are nonzero, 


- - Q 1^1 _ 

cy j'j' ^ j' I -L 

V dr , 


<7ee — Si 


r — R 


(29) 

(30) 


Here, Sr and Sg are the moduli of radial and circumferential hbers, respec¬ 
tively. At steady state, the inertial term in the momentum equations of huid 
and structure vanishes, and we have V-CTe = 0. Therefore, along the r-direction 
we get 


da, 


e,rr ^ ^e,rr ^e,66 


dr r 

This gives us the equation for pe as 


diyGrr Pe) ^ ^rr ^08 _ g 

dr r 


(31) 


dr 


darr Grr ^80 _ ^ | 1 

~dfr~ ^ r “ I r 


IdP 
r dr 




(32) 


We impose zero-traction boundary condition on the four lateral surfaces of 
the huid domain. At steady state, this implies that the exterior huid pressure 
is zero, so that traction continuity at the exterior huid-solid interface implies. 


0 = n • (-Pel -I- cr) • n = -Pe(ro) + Cr^r(?'o)- (33) 

The solution to eq. (32) for pe with the boundary condition given by eq. (33) 
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reads as 


Pe{x) = Sr 


1 



IdR 
r dr 



The inner pressure, Pinner is obtained from eq. (27) as 


Pinner Pe(^i) S, 

dR 


= Sr 


dr 


in) 


dR 


dR. 


dr 


+ / Sr 

Jto 


1 

r 


1 dR 
r dr 




dr. 


(35) 


Under dilation and in the absence of shearing motions, the radial fibers will 
become compressed. For compression-resistant fibers such as those used in this 
model, this is an unstable configuration, and the minimum energy configura¬ 
tion is obtained when the layers of the tube rotate to release the energy of 
compression. In this configuration, the radial fibers do not contribute to the 
elastic stress. Therefore, by letting Se = S and Sr = 0, the fluid pressure at 
the inner surface r* becomes 


p. = 

inner 




-s(--^) dr. 


Substituting eq. (28) into eq. (36), we obtain the explicit form of Pir 


(36) 

as 


p. — Q 
^ inner ^ 




log(ri + Jrf - C) 


\7Uc 


(37) 


log {r^ + Jrl- C] + 


rl-C 


- log ( ^ 

V I n 


Here C = rf — R?, ri = Ri + U{ri), U{ri) denotes the radial displacement of 
the inner surface of the tube, and Tq = JR^ — Pf -|- rf. 
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Fig. 1. (a) Schematic of patches (hatched areas) associated with radial springs (blue 
lines) and circumferential springs (green curves) in (r, 9) plane. The red curves de¬ 
note the boundaries of the tube: the inner (r = a) and the outer surface {r = b). In 
(r, 9) plane, the patch of a radial spring connected by nodes (r, 0, z) and (r +Ar, 0, z) 
is (r, r + Ar) x (9 — O.5A0, 9 + O.5A0). The patch of a circumferential spring con¬ 
nected by the nodes {r^9^z) and (r, 0 -h A9^z) is (r — A^,r -h A^) x {9^9 -\- A9) 
for a < r < 6; (r — A^,r) x {9,9 + A9) for r = 6 and (r, r -h A^) x {9,9 + A9) 
for r = a. Here Aj. = r — ^r(r — Ar) and A^ = —r -h y/r(r -h Ar) are such that 
the circumferential patch is partitioned into two equal volumes, (b) Schematic of 
patches (hatched areas) associated with radial springs (blue lines) and axial springs 
(magenta lines) in (r, z) plane. The red curves denote the boundaries of the tube: 
the inner (r = a), outer {r = b), lower {z = 0) and the upper surface {z = 1). In (r, z) 
plane, the patch of a radial spring connected by nodes (r, 9, z) and (r -h Ar, 9, z) 
is (r, r -h Ar) x {z — 0.5Az, z -h O.SAz) for 0 < z < /; (r, r -h Ar) x {z — 0.5Az, z) 
for z = I and (r, r -h Ar) x {z,z + 0.5Az) for z = 0. The patch of an axial spring 
connected by nodes (r, 9, z) and (r, 9,z + Az) is (r — Aj, r -h A^) x {z,z + Az) for 
a < r < b] {r — Al,r) X {z, z + Az) for r = 6 and (r, r -h A^) x {z,z + Az) for r = a. 
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(a) (b) 

Fig. 2. (a) Schematic of patches (hatched areas) associated with axial beams in 
(r, z) plane. The red lines denote the boundaries of the tube: the top z = I and the 
bottom z = 0. The blue patch is associated with a beam that has a node on the top 
or bottom surface of the tube and its volume is: V = l.briAzAOAr for a < r < 6; or 
V = O.75riAzA0Ar^ if (r — a){r — h) — 0. The magenta patches are associated with 
a beam which does not have a node at the top or bottom surface of the tube. The 
associated patch volume is: V = ViAzAOAr for a < r < 6; or 1/ = O.Sr^AzA^Ar, if 
(r — a){r — h) — 0. (b) Schematic of patches (hatched areas) associated with radial 
springs (blue lines) and circumferential springs (green curves) in (r, 9) plane with 
nonuniform A9. The red curve denotes the interface between the two layers with 
different spacing in 9. Compared with the case of uniform A0, only circumferential 
springs located on the interface need special treatment. 



Longitudinal Muscle 
Vagus Nerve 

Circular Muscle 
Auerbach’s Plexus 
Meissner’s Plexus 


Fig. 3. Schematic of the multiple layers of the esophagus. The inner most layer is 
the mucosal layer (including mucosa and sub mucosa), which is highly folded at rest; 
the outer layers are muscle layers including circular muscle layer and longitudinal 
muscle layer (reproduced with permission from Kahrilas [26]). 
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Fig. 4. Schematic (not drawn to scale) of the computational domain consisting of 
the elastic esophagus (red) and a viscous fluid (green). The elastic esophagus, a 
cylindrical tube with its top end fixed, is immersed in the background fluid in our 
3D computational model. The upper esophagus is initially filled with a bolus, and 
the lower part is filled with a thin liquid layer in the lumen. Traction-free boundary 
conditions are applied to all surfaces of the rectangular computational domain. 



Fig. 5. Axial velocity in the plane ^ = 0 at different times for Case 1 in Section 5.1. 
Only the inner mucosal (IM) layer (white) of the esophagus is shown to better 
visualize the inside bolus. 
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x(mm) t = 0.6s t = 1.2s t = 1.8s t = 2.4s 
t = Os 


Fig. 6. Pressure in the plane ^ = 0 at different times for Case 1 in Section 5.1. Only 
the inner mucosal (IM) layer (white) of the esophagus is shown to better visualize 
the inside bolus. 
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Fig. 7. Kinematics of esophageal layers at four different stages: at rest (t=0 
s); at dilation (t=1.14 s); at contraction (t=1.28 s); and at relaxation 
(t=2.4 s) for Case 1 in Section 5.1. Purple, blue, magenta, green and or¬ 
ange meshes from the inside to the outside, denote the IM, OM, IF, CM 
and LM layers, respectively. (Upper) Side view of a section of the esophagus 
within the box: (—7 mm, 7 mm) x (—0.2 mm, 0.2 mm) x (75 mm, 175 mm); 
(Lower) top view of a section of the esophagus within the box: 
(—7 mm, 7 mm) x (—7 mm, 7 mm) x (119.5 mm, 120.5 mm). Because of 
complex kinematics of the esophageal structure, the apparently overlapping or 
missing springs in the above figures are actually a consequence of the out-of-plane 
motions of the springs. 
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Fig. 8. The cross-sectional area (CSA) of the bolus and the esophageal components, 
and the lumen pressure along its central line: x = 0, 2 / = 0, at t = 1.2 s for Case 1 
in Section 5.1. 
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^t^Os ^ = t = 1.2s t = 1.8s t = 2.4s 


Fig. 9. Axial velocity in the plane ^ = 0 at different times for the Case 2 in Sec¬ 
tion 5.2. Only the inner mucosal (IM) layer (white) of the esophagus is shown to 
better visualize the inside bolus. 
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Fig. 10. Pressure in the plane ^ = 0 at different times for the Case 2 in Section 5.2. 
Only the inner mucosal (IM) layer (white) of the esophagus is shown to better 
visualize the inside bolus. 
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Fig. 11. Kinematic information of esophageal layers at four different stages: 
at rest (t=0 s); at dilation (t=1.14 s); at contraction (t=1.48 s); at relaxation 
(t=2.4 s), for the Case 2 in Section 5.2. Purple, blue, magenta, green and 
orange meshes from the inside to the outside, denote the IM, OM, IF, CM 
and LM layers, respectively. (Upper) Side view of a section of the esophagus 
within the box: (—7 mm, 7 mm) x (—0.2 mm, 0.2 mm) x (75 mm, 175 mm); 
(Lower) top view of a section of the esophagus within the box: 
(—7 mm, 7 mm) x (—7 mm, 7 mm) x (119.5 mm, 120.5 mm). 
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Fig. 12. The cross-sectional area (CSA) of the bolus and the esophageal components, 
and the lumen pressure along its central line: x = 0, ^ = 0, at t = 1.2 s for Case 2 
in Section 5.2. 
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(a) 


(b) 


Fig. 13. (a) Schematic of two families of helical fibers (the blue and red lines, re¬ 
spectively) and the patch (the hatched area) in (z, 6) surface with radial coordinate 
r. The patch volume associated with the helical spring of each family is rAOArAz^ 
if a < r < 6; or O.brAOArAz^ if r = a or r = 6, where a and b are the r-coordinate 
of the inner surface and outer surface, respectively, (b) Illustration of helical fibers 
running on the inner-most surface of the mucosal layer. 



x (mm) t = 0.6s t = 1.2s t = 1.8s t = 2.4s 
t = Os 


Fig. 14. Axial velocity in the plane, ^ = 0 at different times for the Case 3 in 
Section 5.3. Only the inner mucosal (IM) layer (white) of the esophagus is shown 
to better visualize the inside bolus. 
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x(mm) t = 0.6s t = 1.2s t = 1.8s t = 2.4s 


Fig. 15. Pressure field in the plane, y = 0 at different times for the Case 3 in 
Section 5.3. Only the inner mucosal (IM) layer (white) of the esophagus is shown 
to better visualize the inside bolus. 



Fig. 16. Kinematic information of esophageal layers at four different stages: 
at rest (t=0 s); at dilation (t=1.14 s); at contraction (t=1.36 s); at relaxation 
(t=2.4 s), for the Case 3 in Section 5.3. Purple, blue, magenta, green and 
orange meshes from the inside to the outside, denote the IM, OM, IF, CM 
and LM layers, respectively. (Upper) Side view of a section of the esophagus 
within the box: (—7 mm, 7 mm) x (—0.2 mm, 0.2 mm) x (75 mm, 175 mm); 
(Lower) top view of a section of the esophagus within the box: 
(—7 mm, 7 mm) x (—7 mm, 7 mm) x (119.5 mm, 120.5 mm). 
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Fig. 17. The cross-sectional area (CSA) of the bolus and the esophageal components, 
and the lumen pressure along its central line: x = 0, ^ = 0, at t = 1.2 s for Case 3 
in Section 5.3. 
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